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1.  Executive  Summary 


This  report  documents  a  summary  of  the  work  accomplished  during  the  first 
six  months  under  the  project  entitled  “Identification  of  Airborne  Radar  Targets” 
(Contract  #  N00014-86-K-0202).  During  this  period  of  the  renewed  contract,  the 
progress  toward  realizing  the  goals  set  forth  in  the  statement  of  work  has  been  rapid, 
and  has  produced  a  number  of  key  interim  results.  In  particular,  we  have: 

1.  Developed  techniques  to  estimate  and  parameterize  canonical  scattering  cen¬ 
ters  of  a  radar  object.  The  results  of  this  technique  are  compared  with  those 
obtained  using  an  ARMA  estimation  technique  for  single  polarization  data. 

2.  Developed  and  evaluated  the  performance  of  a  signal  processing  technique  that 
estimates  the  polarization  characteristics  of  target  scattering  centers.  This 
technique  incorporates  polarization-diverse  data  to  obtain  a  super-resolved 
down-range  cross-section  of  the  target  showing  polarization  ellipses  of  the 
predominant  scattering  centers. 

3.  Evaluated  the  feasibility  of  using  bispectrum  processing  techniques  for  the 
purpose  of  radar  target  characterization  and  identification.  The  results  of 
this  study  have  identified  advantages  of  bispectrum  processing  for  frequency- 
-di-persive  scatterers  and  indicate  the  possibility  of  identifying  targets  from 
the  estimated  bispectral  responses. 

4.  Developed  improved  techniques  for  estimating  the  statistical  parameters  of 
observed  symbol  strings  for  the  purpose  of  target  identification  based  on  struc¬ 
tural  pattern  recognition.  In  addition,  we  have  developed  a  iormuiation  of  the 
target  identification  algorithm  that  resembles  a  tree  classifier.  Reformulating 
the  algorithm  in  this  way  reveals  the  structural  nature  of  the  classifier  to  the 
extent  that  a  relative  evaluation  of  the  discrimination  capability  of  particular 
symbols  may  be  possible. 

The  remainder  of  this  report  is  organized  into  chapters  that  more  fully  discuss 
progress  and  ongoing  research  in  the  areas  mentioned  above. 
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2.  High  Resolution  Radar  Target  Modeling  Using  ARMA  Models 

This  chapter  discusses  research  high  resolution  modeling  of  radar  data  to  ob¬ 
tain  target  impulse  responses.  The  techniques  are  parametric,  and  rely  on  scatter¬ 
ing  center  location  estimation,  and  parameterizing  the  scattering  center  by  various 
attributes.  This  chapter  focuses  on  single  polarization  data.  Two  methods  are 
considered.  One  is  a  canonical  scattering  center  analysis,  which  employs  nonlin¬ 
ear  minimization  techniques  to  obtain  the  estimate.  The  other  approximates  this 
canonical  scattering  model  by  an  Autoregressive  Moving  Average  (ARMA)  model, 
and  uses  ARMA  estimation  techniques  derived  from  the  signal  processing  literature. 
The  latter  method  is  more  computationally  efficient  than  the  former  method. 

We  begin  with  a  discussion  of  canonical  scattering  centers.  We  next  present  the 
ARMA  modeling  method,  and  briefly  describe  the  algorithm  which  is  used.  Next 
we  present  simulation  results  obtained  using  compact  range  data  measurements 
of  model  aircraft.  Finally,  we  discuss  the  canonical  scattering  center  estimation 
method. 

2.1  Canonical  Scatterers 

It  is  well-known  [1]  that  scattering  of  electromagnetic  waves  at  high  frequencies, 
i.e.  at  frequencies  such  that  the  wavelength  is  small  compared  to  the  scatterer,  can 
be  modeled  as  a  set  of  independent  scattering  centers.  There  are  several  varieties 
of  canonical  scattering  centers  used  in  this  process;  typical  examples  are  the  point, 
edge  and  plate  scatterers  and  the  corner  reflector  [1]. 

If  a  target  consisting  of  m  scatterers  is  illuminated  using  a  CW  radar,  then  the 
backscattered  field  received  by  a  linearly  polarized  antenna  will  give  a  signal  of  the 
form  [2]: 

»*  =  (2.1, 

■=1 

where  A',  is  the  amplitude,  r,  the  relative  range,  and  f,  the  type  of  the  scatterer; 
fk  is  the  kfi1  frequency,  and  m  is  the  number  of  scatterers.  For  the  canonical 
scatterers  in  [2],  tx  takes  on  the  following  values 
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/»€{—•  -0.5,0, 0.5,1}  (2.2) 

In  addition,  for  stepped  frequencies,  we  can  represent  fk  by 

fk  =  fo  +  Mj,  k  =  0,1,2,. . .  ,N  —  1  (2.3) 

where  /0,  6 j  and  N  are  given  parameters  of  the  system  under  consideration.  The 
frequency  /0  represents  the  lowest  frequency  the  radar  system  transmits  and  8j 
is  the  size  of  the  frequency  step.  Estimation  of  the  parameters  in  (2.1)  from  N 
measurements  leads  to  a  canonical  scattering  description  of  the  target;  see  [2], 

This  model,  while  based  on  the  physical  properties  of  the  scattering  process, 
has  the  disadvantage  of  being  non-linear.  A  non-linear  estimator  typically  requires 
more  data  to  converge  to  within  a  given  distance  of  the  correct  estimate,  and  also 
suffers  from  the  possibility  of  convergence  to  a  false  (local,  but  not  global)  minimum, 
leading  to  incorrect  estimates. 

We  are  currently  investigating  two  estimators  for  this  model.  The  first  estimator 
uses  the  model  (2.1)  directly,  and  estimates  the  parameters  of  that  model  through 
a  multi-step,  non-iterative  algorithm.  If  desired,  the  quality  of  the  resulting  esti¬ 
mates  can  be  improved  further  by  using  the  multi-step  estimate  as  an  initial  guess 
for  an  iterative  non-linear  least-squares  estimator.  The  second  estimator  uses  an 
approximate  model  of  (2.1)  that  allows  linear  estimation  techniques,  specifically  an 
ARM  A  algorithm.  The  second  algorithm  is  more  computationally  efficient,  and  has 
been  studied  more  extensively. 

The  approximate  model  for  linear  estimation  is  obtained  by  replacing  the  (j2nfkY 
factor  in  (2.1)  with  pk .  By  suitably  adjusting  the  amplitude,  the  approximate  model 
can  be  made  to  match  the  original  model  if  the  bandwidth  involved  is  not  too  large. 
Thus  we  arrive  at: 

m 

V*  =  T.d,r1  (2.4) 

1=1 

where  d ,  is  in  general  not  equal  to  K,  in  (2.1)  because  of  the  approximation  involved. 
The  variable  p,  can  be  found  from: 

p,  =  (2.5) 
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In  practice  the  absolute  range  is  not  used;  instead,  r;  represents  the  range  relative  to 
a  zero  range  reference  point.  R  then  represents  the  maximum  unambiguous  range 


It  should  be  noted  that  for  t  =  0  the  correspondence  is  exact  (p  will  be  1). 

The  energy  associated  with  a  single  scattering  center  is  the  square  of  the  absolute 
value  of  its  contribution  to  the  amplitude,  or 
N-i  i  _  Q2N 

p,  =  z  i<y2ip.r = (2-v 

,=o  1  -  Pi 

Equations  (2.4)  and  (2.6)  will  be  used  to  form  the  approximate  parametric  model 
of  target  scattering. 


2.2  ARMA  Models 


In  this  section  we  will  discuss  estimation  techniques  for  the  approximate  model. 
For  a  more  detailed  analysis  of  this  model  and  the  corresponding  estimation  tech¬ 
niques,  see  [3]  and  [4]. 

Equation  (2.4)  can  be  formed  as  the  solution  of  a  linear  difference  equation, 


m 

yk  +  J2a>yk+' 

m-1 

•=  ^  1  Ci6k+i 

(2.7) 

i=i 

i=0 

A(q~1yk 

=  C{q-')6k 

(2.8) 

where  <5*.  is  the  unit  impulse  sequence.  It  can  be  shown  [5]  that  the  {p,}  in  equation 
(2.4)  are  the  roots  of  the  A(z)  polynomial,  and  the  coefficients  {d,}  the  partial 
fraction  residues  of  C(z)  and  A(z).  It  should  be  observed  that  this  means  that  the 
radar  cross  section  as  a  function  of  relative  range  can  be  found  by  taking  the  square 
of  the  absolute  value  of 


r/R) 

^r-  =  yl(ej>(l-2r/ft) 

f*  dj 

2-s  ej*(l-2r/R )  _  p.' 


0  <  r  <  R 
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By  observing  that  for  k  >  0  the  right  hand  side  in  (2.7)  equals  zero,  we  ran  find 
a  linear  least  squares  estimator  for  the  A(z)  polynomial: 


yM+ 1 

Vm  ■ 

y\ 

O'M 

VM+ 2 

y.vr+ 1  • 

V2 

'■ 

ay 

y.v-i 

VN-2  ■ 

■  yN-M-i 

1 

where  M  is  some  integer  at  least  equal  to  m  and  where  equality  will  hold  for  noiseless 
data  that  fits  the  model  perfectly.  Appropriate  ways  for  choosing  M  are  discussed 
in  [4]. 

The  performance  of  this  estimator  can  be  improved  through  the  use  of  the 
singular  value  decomposition  (SVD)  as  follows:  find  the  SVD  of  the  matrix  Y  on 
the  left  hand  side  of  (2.9);  let 


O  i  >  O 2  >  .  .  .  >  <7 M+ 1  ^  0. 


be  the  resulting  singular  values.  Set  ern+i . .  .  a M  +  1  equal  to  zero,  reconstruct  1' 
and  solve  (2.9)  with  the  new  Y.  The  value  of  n  is  chosen  to  be  such  that  the 
first  large  drop  in  the  value  of  the  {oq}  falls  between  an  and  <7n+j.  This  choice 
is  motivated  by  the  ability  of  the  SVD  to  partially  reject  the  effects  of  noise  and 
unmodeled  dynamics  on  the  estimate.  This  issue  is  discussed  in  further  detail  in 


[4]- 

Once  the  A(z)  polynomial  has  been  found,  the  {p*}  can  be  estimated  through 
standard  polynomial  rootfinding  techniques.  Given  the  {pj,  the  {d,}  can  be  found 
from 


■  p?  • 

-  go 

Vo 

.  p?" 1 

..  jj.jg 

<R, 

dN_  i 

.  Vn-i  . 

(2.10) 


where  p,  i«  the  e  '•’mate  for  p,  found  above. 

Above  we  -".puted  an  estimate  of  a  polynomial  A\f(z)  of  order  M  >  m,  and  a 
corresponding  set  of  M  zeros  {p,}^.  We  also  computed  amplitudes  corresponding 
to  these  zeros.  The  next  step  is  to  separate  the  m  zeros  that  correspond  to  the 
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radar  target  (signal  zeros)  from  the  remaining  spurious  zeros.  The  zero  separation 
process  relies  on  the  observation  that  since  the  spurious  zeros  correspond  to  modes 
that  are  not  present  in  the  data,  their  estimated  energy  should  be  low.  To  take 
advantage  of  this  property,  we  use  the  following  strategy.  First,  the  amplitudes 
{ c/x } ! ,  corresponding  to  all  M  zeros  are  found  using  equation  (2.10).  Then,  the 
energy  of  each  mode  is  computed  using  equation  (2.6).  The  m  zeros  whose  cor¬ 
responding  energy  is  the  highest  are  kept  and  the  remaining  ones  are  discarded. 
Then,  the  energy  estimates  for  the  m  kept  zeros  can  be  refined  by  recomputing  the 
amplitude  estimation  using  (2.10)  with  only  the  m  selected  zeros.  However,  since 
the  energies  for  the  discarded  zeros  is  low,  the  difference  between  the  two  estimates 
should  generally  be  small,  so  this  refinement  step  may  be  omitted. 

2.3  Experimental  Results 

In  this  section  we  describe  results  obtained  with  both  types  of  algorithm.  We 
first  desciibe  results  obtained  by  applying  the  ARM  A  algorithm  to  data  obtained 
from  the  ESL  Compact  Range  and  then  present  results  from  the  application  of  the 
canonical  estimator  to  simulated  data. 

2.3.1  Results  from  Applying  ARMA  Techniques  to  Compact  Range 
Data 

The  main  aim  of  this  section  is  to  demonstrate  that  the  ARMA  type  algorithm 
presented  above  is  robust  to  changes  in  model  order  M  and  number  of  singular 
values  kept  n,  performs  well  in  the  presence  of  noise  and  achieves  consistent  results 
across  different  frequency  ranges.  On  the  other  hand,  the  algorithm  readily  detects 
differences  between  aircraft  and  is  able  handle  changes  in  aspect  angle. 

Figure  2.1  shows  estimated  responses  for  a  1:130  model  of  the  Concord  using 
A'  =  30  data  points  taken  evenly  spaced  from  18-26  GHz  (corresponds  to  138 
200  MHz  on  a  full-size  Concord)  and  model  orders  M  =  5  and  10,  keeping  n  =  3 
singular  values  and  m  =  2  zeros. 

Figure  2.2  shows  estimates  for  M  =  10  order  model  from  the  same  data  as 
above,  keeping  n  =  3,5  and  7  singular  values  and  m  =  2  zeros. 
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Figure  2.3  shows  estimates  for  M  =  7,  n  =  5  and  m  =  2  but  with  Gaussian 
noise  added  at  10  dB  SNR. 

The  above  figures  demonstrate  by  their  similarity  that  the  algorithm  is  robust 
with  respect  to  the  effects  of  noise,  choice  of  order  and  choice  of  the  number  of 
singular  values  to  keep. 

Figure  2.4  shows  estimated  responses  for  the  Concord  taken  from  N  =  20  data 
points  in  a  6.25-12  GHz  range  (48-92  MHz  for  a  full-scale  target)  using  M  =  10, 
n  =  7  and  m  =  4.  We  have  used  two  data  sets  here,  one  taken  at  an  aspect  angle  of 
zero  degrees,  the  other  at  20  degrees.  The  figure  clearly  shows  the  relation  between 
the  two  estimated  responses.  We  can  see  that  the  second  large  peak  has  split  into 
two  peaks  at  20  degrees,  and  all  the  scattering  centers  are  a  bit  closer  to  each  other, 
as  a  result  of  the  apparent  foreshortening  of  the  target  in  the  downrange  direction. 

For  reference,  Figure  2.5  gives  the  result  of  an  IFFT  performed  on  the  18-26  GHz 
data,  and  Figure  2.6  gives  a  properly  scaled  picture  of  the  Concord. 

2.3.2  Direct  Estimation  of  the  Canonical  Scatterer  Model 

In  this  section  we  detail  an  algorithm  that  directly  estimates  the  {A';},  {n}  and 
{t{}  parameters  of  equation  (2.1).  The  order  of  estimation  is  first  the  {r,},  then  the 
{t, }  and  finally  the  { K, }  parameters. 

The  {r,}  parameters  are  estimated  in  the  same  way  as  as  for  the  approximate 
model  in  the  previous  section.  We  then  find  the  { ti }  parameters  by  applying  the 
following  three  step  procedure  for  all  scatterers. 

1.  For  the  scatterer,  modulate  the  sequence  {y^}  with  e  *=  fk ,  thus  trans¬ 
lating  the  scatterer  to  relative  range  zero. 

2.  Apply  a  low-pass  filter  to  the  modulated  data  y*.  Note  that  this  is  a  low-range 
filter,  not  a  low-frequency  filter.  This  filtering  operation  should  leave  the 
scatterer  as  the  only  component  of  the  filtered  data.  Note  that  it  is  important 
that  the  filter  have  linear  phase  for  the  third  step  to  succeed. 

3.  The  filtered  data  can  now  be  described  by 
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Figure  2.5:  Concord,  1 
degree  aspect  angle. 
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Figure  2.6:  Scaled  Picture  of  the  Concord. 
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From  which  it  follows  that 

Ln  y'k  =  LnKi  +  t.Ln  (j2nfk)  (2.11) 

where  Ln  is  the  main  branch  of  the  complex  natural  logarithm.  The  param¬ 
eters  {f,}  can  be  found  from  a  least  squares  solution  to  equation  (2.11).  and 
finding  the  least  squares  solution. 

Finally,  the  {KJ  can  be  found  using  a  least  squares  method  in  equation  (2.1), 
using  the  estimates  for  {r,}  and  {f,}  found  above. 

If  desired,  the  estimates  found  with  the  aid  of  the  multi-step  procedure  above  can 
be  used  as  an  initial  guess  for  a  non-linear  least  squares  estimator.  The  accuracy  of 
the  estimates  for  simulated  data  has  been  sufficient  that  this  has  not  been  necessary 
so  far,  but  we  anticipate  that  such  refinement  may  be  needed  for  real  data. 
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3.  Exponential  Modeling  of  Fully  Polarized  Radar  Return 


The  horizontal  and  vertical  radar  return  from  a  target  illuminated  with  stepped 
frequency  circularly  polarized  radiation  can  each  be  modeled  a  sum  of  complex 
exponentials.  The  scattering  points  of  the  target  are  thus  the  poles  of  these  models, 
and  amplitudes  associated  with  each  pole  represent  the  scatterer’s  horizontal  and 
vertical  return.  One  approach  is  to  use  this  modeling  on  both  the  horizontal  data 
and  the  vertical  data  independently.  This  model  can  then  be  used  to  identify  the 
target. 

However,  it  is  reasonable  to  assume  that  if  a  scattering  center  appears  in  the 
horizontal  return  then  it  is  also  present  in  the  vertical  return.  Thus  when  estimat¬ 
ing  the  poles  of  the  model,  both  the  horizontal  and  vertical  information  are  used 
simultaneously  to  estimate  the  ranges  of  the  scattering  centers.  To  make  the  esti¬ 
mate  more  resilient  to  noise  we  employ  Singular  Value  Decomposition  (SVD)  on  the 
associated  linear  prediction  matrix.  Also,  since  scatterers  can  be  expected  to  cause 
peaky  responses,  poles  in  the  linear  prediction  model  which  are  not  within  a  certain 
distance  from  the  unit  circle  are  eliminated.  Finally,  we  estimate  the  horizontal  and 
vertical  amplitude  characteristics  independently  using  the  appropriate  polarization 
data. 

Once  the  model  has  been  determined,  the  angle  at  which  the  poles  appear  in  the 
z-plane  represent  the  downrange  location  of  each  of  the  scatterers.  The  horizontal 
and  vertical  amplitudes  associated  with  each  pole  contain  the  information  about 
that  scatterer’s  polarization  effect  on  the  incident  radar  wave.  These  polarizations 
are  in  general  elliptical  and,  along  with  the  scatterer’s  range,  can  be  used  to  identify 
the  target. 

3.1  Modeling  the  Scatterers 

Assume  that  two  sets  of  N  data  points  are  generated  from  M  complex  exponen¬ 
tials  with  white  Gaussian  noise  added,  where  each  set  represents  the  horizontal  and 
vertical  return  from  a  left  circularly  polarized  incident  wave,  and  svi  respectively. 
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Shi{n) 

\1 

=  £ 

O-hk 

pnk  + 

wh(n) 
wv(n)  _ 

_  svi(n)  _ 

k=  l 

-a 

n  =  1, . .  ,  N. 


(3.1) 


Where  pk  is  the  kth  common  pole,  a^k  is  the  horizontal  amplitude  associated  with 
that  pole,  and  avk  is  the  vertical  amplitude  associated  with  the  same  pole.  The 
backward  linear  prediction  equations  can  then  be  set  up  as  follows[6][7][8j. 


Shi(2) 

su(  3)  •  •  ■  sm(L  +  1) 
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Shi{  3) 

^/i/(V  —  L  +  1) 

.  su(N) 

'  k  ' 

shi(N  -  L ) 

Svt(  2) 

Svi{  3) 

Svi(  3)  svi(L  +  1) 

.  h  . 

Mi) 

.  svl(N  -L  +  1) 

.  Sv,(N) 

.  svl(N  -  L) 

(3.2) 


or 


Sb=  -s 


(3.3) 


where  L  is  the  order  of  prediction,  and  b  is  the  coefficient  vector  of  the  polynomial 
B(z)  given  by 


B(z)  —  1  +  b\Z  +  . . .  +  biz 


(3.4) 


Note  that  both  sets  of  data  are  used  simultaneously  to  estimate  a  single  set  of 
prediction  coefficients. 

The  solution  involves  forming  the  matrix  [s  :  5],  performing  a  Singular  Value 
Decomposition  (SVD)  on  it,  then  truncating  all  but  the  first  M  singular  values  to 
arrive  at  a  better  estimate  fs  :  S'] .  Next,  b  is  found  using  the  pseudoinverse  of  S, 


(3.5) 


b  =  -5+5. 

Now  the  estimated  poles  can  be  determined  as: 

1 


P.  = 


root 


.  (£(*)) 


i  —  1 ,  •  ■  •  i  L 


(3-6) 
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Since  scatterers  result  in  peaky  responses,  those  poles  which  do  not  lie  within  a 
given  annular  region  about  the  unit  circle  can  be  eliminated  from  the  model.  The 
following  criterion  has  been  found  to  work  well  for  radar  data: 


1 


I5o  <  *1"  <  100 


(3.7) 


We  keep  only  those  poles  in  (3.6)  which  satisfy  (3.7).  Once  these  V  poles  have  been 
determined,  the  amplitude  equations  for  both  the  horizontal  and  vertical  compo¬ 
nents  can  be  formed. 

*^/l/(l)  £„/(!) 


Pi 


Pu 


or 


L  Pi 


PA  =  Sa 


Pu 


d/ii 

av  l 

.  <*KL 

Q-vL 

Shi(N)  svi(N ) 


(3.8) 


(3.9) 


Total  least  squares  can  now  be  used  to  estimate  the  amplitudes. 

A  =  (p"pylPHsa  (3,io) 

3.2  Determining  Scatterer  Locations  and  Associated  Polarizations 

Now  that  we  have  ahi,  dvl,  and  pt,  for  each  scatterer  i,  we  can  determine  its 
range  r,  and  its  polarization  in  terms  of  tilt  fy,  ellipticity  e,,  and  major  axis  .4,. 
They  are  given  by  the  following: 


II 

(3.11) 

i  i 

r,  =  -tan 

(tan  (2-y.)  cos  (d,)) 

(3.12) 

-  1  -i  1 

t.  —  -sin 

O  1 

{ sin  (27.)  \ 

fc\  J 

(3.13) 
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V  sm  (<5,J  ) 
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sin  (St)  |sec (f,)| 

V  / 

Ai--rJLL 

(3.14) 

_I _ 2tan(t,)cos{6,)  tan^(r,) 

El  EuBh,  +  El 
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where 


7*  =  tan  1 

(fd 

(3.15) 

\EuJ 

—  Lat 

Lahx 

(3.16) 

E\t  — ■ 

(3.17) 

E'2i  —  |^vi| 

(3.18) 

The  ranges  of  a  target’s  scatterers  and  their  respective  polarization  characteris¬ 
tics  can  thus  be  used  to  identify  the  target. 


3.3  Simulation  Results 


The  above  estimation  procedure  was  applied  to  data  obtained  from  compact 
range  measurements  of  several  small  aircraft  models;  in  particular,  we  used  scale 
models  of  the  707,  727,  747,  DC10,  and  the  Concord.  The  ESL  compact  range 
was  used  to  measure  horizontal  transmit  with  horizontal  and  vertical  receive  and 
vertical  transmit  with  horizontal  and  vertical  receive,  Shh(f),  $vh(f)  =  Shv(f),  and 
svv(f)  respectively.  The  measurements  were  made  for  frequencies  /  between  1.5  and 
12  GHZ  at  50  MHz  steps,  thus  resulting  in  211  data  points  for  each  polarization. 
From  this  data,  Shiif )  and  svi(f )  are  found  using  the  following  transformation: 


su(f)  ' 

Shh(f)  Shv(f) 

’  1  ’ 

_ » 

.  Svh(f)  Svv{f)  _ 

.  j  _ 

_1_ 

V5 


(3.19) 


The  following  simulations  were  done  with  N  =  80  data  points  from  the  upper 
half  of  frequency  range.  In  all  cases,  the  nose  on  aspect  angle  data  was  used. 

The  models  range  from  about  23  to  43  cm  in  length,  so  an  unambiguous  range 
R  =  75 cm  is  needed  to  get  resolution  between  the  model’s  individual  scatterers. 
Since  Af  =  2^,  where  c  is  the  speed  of  light,  for  this  R,  Af  =  200Mhz.  This 
corresponds  to  every  fourth  point  in  the  data  set.  To  allow  use  of  all  80  points, 
the  backward  prediction  equations  were  written  so  as  to  interleave  four  sets  of  20 
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points  each.  The  amplitudes  associated  with  each  pole  similarly  used  the  four  sets 
of  data. 

Finally,  model  order  and  the  number  of  singular  values  to  keep  were  chosen  as: 
L  =  10  and  M  =  5.  The  results  for  each  plane  with  no  noise  added  to  the  data 
appear  in  Figs.  3. 1-3.5.  Each  scatterer  is  represented  by  an  ellipse  which  describes 
its  polarization  effect  at  the  appropriate  range.  The  scale  is  in  meters  and  is  for 
ellipse  location.  The  ellipse  factor  of  .2  can  be  used  to  multiply  the  range  scale  to 
see  ellipse  magnitudes.  A  scale  plan  view  for  each  plane  also  appears  in  order  to 
match  up  ellipses  with  their  physical  scatterer. 
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Noise  was  then  added  for  SNRs  of  10 dB,  5 dB,  and  0 dB  to  the  DC10  data.  Five 
simulations  with  different  noise  realizations  for  each  SNR  appeax  in  Figures  3.6- 
3.8.  Since  the  noise  tends  to  make  discrimination  of  ellipses  difficult,  a  plot  of  range 
vs.  the  length  of  the  hypotenuse  between  each  ellipse’s  major  and  minor  axis  also 
appears  in  each  figure. 
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The  simulation  results  demonstrate  that  the  procedure  produces  target  signa¬ 
tures  which  are  unique  and  identifiable  and  that  it  stands  up  well  in  noise,  par¬ 
ticularly  in  finding  the  existence  and  range  of  a  scatterer.  This  is  due  to  the  fact 
that  both  sets  of  data,  s^i  and  svi,  are  used  to  determine  the  poles  of  the  model. 
Even  at  OdB  SNR  the  scattering  centers  are  quite  distinguishable,  particularly  in 
the  hypotenuse  plots.  The  polarization  ellipse  for  each  scatterer,  however,  is  not  as 
robust  since  its  components  are  each  derived  from  half  of  the  data  separately.  Still, 
the  scatterer  center  location  makes  identification  of  the  target  possible  in  these  high 
noise  situations.  The  scatterers  also  appear  at  downrange  locations  which  corre¬ 
spond  with  features  of  the  planes  such  as  wings,  engine  inlets,  rudders,  etc.,  which 
gives  physical  significance  to  the  model  and  helps  to  justify  its  use. 
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4.  A  Study  of  The  Bispectral  Estimates  of  Radar  Signals  And 
Application  To  Radar  Target  Identification 


4.1  Introduction 

This  Chapter  is  concerned  with  the  application  of  bispectrum  estimation  to 
the  analysis  of  radar  signatures.  The  recent  interest  in  bispectrum  estimation  is 
motivated  by  the  broad  area  of  applications  of  higher  order  spectral  estimates, 
especially  the  bispectrum.  Problems  such  as  phase  decoupling,  signal  nonlinearities 
and  skewness  from  Gaussian-like  features  can  be  solved  and  better  understood  using 
bispectrum  estimation.  The  use  of  higher  order  spectra  is  not  an  alternative  for 
lower  order  spectral  analysis  but  rather  a  complementary  process.  Many  features 
embedded  in  measured  data  may  not  appear  in  lower  order  spectral  analysis  or 
could  be  misinterpreted  or  misread.  High  order  spectra  may  uniquely  identify  such- 
features  and  present  them  in  a  way  which  leaves  less  room  for  misinterpretation. 

Radar  measurement  data  contain  a  wide  spectrum  of  information  that  corre¬ 
sponds  to  various  scattering  mechanisms.  In  this  study  an  attempt  is  made  to  use 
bispectral  analysis  to  identify  such  mechanisms.  This  attempt  is  motivated  by  the 
fact  that  such  mechanisms  are  often  confused  using  first  order  spectral  analysis. 
Moreover,  radar  targets  at  certain  signal  to  noise  ratio  and  certain  conditions  of 
azimuth  and  elevation  may  have  similar  lower  order  spectral  responses  but  different 
higher  order  spectral  responses.  Thus,  the  bispectrum  can  be  used  to  discriminate 
among  targets  that  differ  in  fine  details  only. 

Therefore,  the  two  main  reasons  for  studying  bispectrum  in  radar  signature  anal¬ 
ysis  is  to  recover  complex  scattering  features  and  to  discriminate  among  targets  that 
differ  in  fine  details  rather  than  major  features.  Also,  in  this  study  we  looked  at 
bispectral  responses  for  frequency  dispersive  scatterers  and  found  some  interesting 
advantages  of  the  bispectrum  under  these  conditions.  Moreover,  due  to  the  correla¬ 
tion  involved  in  bispectrum  estimation,  the  behavior  of  bispectral  responses  under 
conditions  of  high  noise  is  relatively  less  sensitive  to  noise  corruption  than  lower 
order  spectral  responses.  Finally  an  attempt  is  made  to  classify  radar  targets  based 
on  the  corresponding  bispectral  responses. 
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4.2  r'roulem  Fom.ulatlon:  Birange  Profile  Theory 

Radar  scattering  from  an  object  is  often  measured  as  amplitude  and  phase  at  a 
set  of  evenly  spaced  frequencies.  This  permits  a  transformation  to  the  time  domain 
and  thus  the  range  domain  by  taking  an  inverse  Fourier  transform  of  the  data.  The 
result  is  called  the  radar  target,  impulse  response  [9].  In  general  there  is  a  one  to  one 
relationship  between  specific  scattering  mechanisms  and  the  time  such  mechanisms 
appear  in  the  impulse  response.  One  of  the  difficulties  that  this  type  of  analysis  has 
is  that  complex  targets  often  have  multiple  interactions.  The  result  is  that  many  of 
the  terms  in  the  impulse  response  are  related  to  the  interactions  rather  than  due  to 
the  result  of  simple  subcomponent  scattering.  Many  of  these  multiple  interaction 
mechanisms  can  be  identified  as  such  by  the  application  of  the  bispectrum  to  the 
radar  scattering  data. 

Also,  this  study  indicates  that  the  individual  target  scattering  mechanisms  are 
more  easily  observed  using  the  bispectrum  under  conditions  of  high  noise  or  target 
dispersive  effects.  It  will  be  shown  that  these  bispectral  responses  constitute  robust 
features  which  may  be  used  in  radar  target  identification  schemes. 

The  radar  scattering  from  a  complex  target  can  be  effectively  modeled  as  a 
sequence  of  discrete  scatterers.  With  this  model  each  element  of  the  set  of  complex 
radar  scattering  data 

Sp(fj)  =  Ap  exp  ( —j2nRrf]/c).  (4.1) 

where 

♦ 

Sp{fj)  =  pf/lscattering  coefficient  at 
jth  frequency 

Ap  =  amplitude  of  pth  scat  tering 
coefficient 

R,p  —  line  of  sight  range  between 

scatterers  and  radar  zero  phase 
(&:  time)  reference  plane. 
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fj  =  jth  radar  carrier  frequency 
c  =  speed  of  light. 

This  can  be  usefully  converted  to  scalar  form  by  simply  taking  the  real  part  of  the 
scattering  value. 

3?{5p(/j)}  =  Apcos{2irRpfj/c)  (4.2) 

Note  that  when  data  are  represented  in  this  form  for  the  scattering  from  a 
particular  single  point  scatterer  at  a  particular  range  from  the  radar  phase  zero 
reference  plane  (which  may  be  located  near  or  on  the  target),  a  sinusoidal  variation 
in  the  data  as  a  function  of  frequency  results.  The  amplitude  of  the  sinusoid  is 
proportional  to  the  scattering  coefficient  of  the  scatterer,  and  the  periodicity  of  the 
sinusoidal  variation  (as  a  function  of  radar  carrier  frequency  )  of  the  real  part  of 
the  the  scattered  signal  is  proportional  to  the  line  of  sight  distance  of  the  scatterer 
to  the  zero  phase  reference  plane.  If  the  scattered  do  not  interact,  then  the  total 
received  signal  will  have  the  form 

<h(/,)  =  '£Apco6(2nRpfj/c)  (4.3) 

V 

Note  that  computation  of  the  spectrum  of  St  will  yield  values  of  Ap  for  each  Rp. 
Note  on  the  other  hand  that  if  the  scatterers  interact  (multiple  interaction),  then 
a  signal  scattered  from  the  latter  part  of  a  particular  target  will  be  influenced  (or 
modulated)  by  interaction  with  components  of  the  target  closer  to  the  radar.  This 
interaction  with  other  parts  of  the  radar  target  prior  to  scattering  from  a  particular 
object  on  the  radar  target  results  in  modulation-like  effects  in  the  range  domain 
profile  w'hich  can  be  uniquely  identified  when  the  bispectral  estimation  process  is 
applied  to  St. 

The  definition  of  the  bispectrum  is  the  two  dimensional  Fourier  transform  of  the 
third  moment  sequence.  In  this  study  the  bispectrum  is  a  birange  profile  obtained 
by  computing  a  two  dimensional  Fourier  transform  of  the  third  order  autocorrelation 
function  which  is  function  of  frequencies  f\,f2  in  the  two  dimensional  frequency 
domain. 
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fls,(/„/2)  =  EiStiMStif,  +  fiWf,  +  /a)} 


(4.4) 


Bs,(RuRi)  =  EEMJ.)  (4-5) 

A  h 

exp(-j2ir/c[Rlf1  +  -R2/2]) 

where  £{.}  stands  for  the  expected  value  and  i?4|(/i,/2)  is  the  third  moment  se¬ 
quence  of  S(t).  The  birange  is  denoted  by  R2)  where  i?i,i?2  belong  to  the 

two  dimensional  down  range  domain. 

The  bispectral  estimates  in  this  study  are  computed  using  the  indirect  classical 
method  which  is  based  on  computing  an  estimate  of  the  third  moment  sequence. 
Windowing  can  also  applied  to  the  estimates  of  the  autocorrelation  values.  In  this 
study  the  optimum  window  (minimum  bispectrum  bias  supremum)  [10,11]  is  used. 

The  above  definition  can  be  used  as  an  indirect  way  to  estimate  the  bispectrum 
but  to  get  some  kind  of  a  physical  insight  of  what  the  bispectrum  really  means  it 
is  better  to  use  the  following  which  is  considered  as  a  direct  way  of  estimating  the 
bispectrum. 

v(r)  =  E  $(/)  exp  (+j2vjj)  (4.6) 

r2)  =  ^(rl)y(r2)?/(r1  +  r2)  (4.7) 

where  y(r)  is  the  impulse  response  and  Bs,(Ri,  R2)  is  the  direct  estimate  of  the 
birange  and  0  is  a  constant  (0  =  (yj)2.  where  f,  is  the  sampling  frequency  of  the 
radar  frequency  data.  Thus,  if  y(r)  is  the  impulse  response  of  the  radar  backscatter 
then  the  birange  f?s,(i2i,  i?2)  is  nonzero  only  when  y(fi),y(r2),  and  y(r1  -f  r2)  are 
not  equal  to  zero.  This  is  mostly  the  case  when  we  have  an  interaction  term  among 
different  scatterers.  Notice  that  b(r l7rq)  =  0y2{ri)y(2ri)  which  indicates  that  there 
exist  scatterers  at  rq  and  2 rq  respectively.  Also,  it  is  important  to  notice  that; 

6(^,0)  =  y(r1)y(0)y*(r1)  (4.8) 
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which  means  that  for  the  bispectrum  to  exist  at  either  time  (range)  axis  we  need  to 
have  a  scatterer  at  r  =  0  then  the  rt  axis  represents  the  impulse  response,  similarly 
for  the  other  axis.  Thus  if  a  scatterer  exists  at  r  =  0  then  the  bispectrum  at  the  r 
axis  represents  the  impulse  response. 

An  expression  for  the  third  moment  sequence  of  complex  data  can  be  obtained 
from  the  direct  method  as  follows; 

R(w,  n)  -  E{Sm(f)S(f  +  m)S(f  +  n)}  (4.9) 

where 

R(m,n)  =  R(n,m)  (4-10) 

Thus  the  only  symmetry  property  for  the  bispectrum  in  this  case  is 

BSt(Ri,R2)  =  BSt{R2,Ri)  (4-11) 

Moreover,  for  complex  data  we  can  use  more  than  one  form  of  bispectrum  as  sug¬ 
gested  by  Brillinger  and  Rosenblatt  [12]  and  that  is  by  placing  the  conjugate  at 
different  combinations  of  S{f).  However,  the  above  definition  is  most  appropriate 
for  interpretation  purposes.  Also,  the  real  part  of  the  radar  data  could  be  used  only 
to  compute  the  bispectrum  with  little  loss  of  information. 

4.3  Classical  Bispectrum  Estimation  Algorithms 

Classical  methods  for  bispectral  estimation  were  used  in  this  study.  High  res¬ 
olution  technique  that  involve  ARMA  modeling  algorithms  may  be  used  also  to 
obtain  bispectral  estimates.  Classical  bispectral  estimators,  often  called  conven¬ 
tional  estimators  are  divided  into  two  categories-  Indirect  methods  that  follow  the 
definition  of  the  bispectrum  as  the  two  dimensional  Fourier  transform  of  the  third 
order  moment  sequence,  and  Direct  methods  similar  to  the  periodogram  estimate 
of  the  spectral  density. 
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4.3.1  Indirect  Method 


Let  {.Y(l),  X(2), .  ,.,X(N)}  be  the  frequency  data  set.  Then  the  estimate  of 
the  third  moment  sequence  is 


R(m ,  n)  =  —  ^2  X(i)X(i  +  m)X(i  +  n)  (4.12) 

i=3l 

where  sx  =  max  (0,  —  m,  —n)  and  s2  =  min  (N  —  1,7V—  1  —  m,N  —  1  —  n).  Then 
the  estimate  of  the  bispectrum  B(Ri,  f?2)  &s  a  function  of  range  R  is  generated  as: 
L  L 

5(ri,r2)  =  Y2  R(m,n)W(m,n)exp(—j4TT/c(r1m  +  r2n)  (4.13) 


where  L  <  N  —  1  and  W(m,  n)  is  a  two  dimensional  window  function.  The  optimum 
window  is  used  in  this  study  and  it  can  be  generated  from  a  one  dimensional  window 
as: 


W(m,  n) 
d(m ) 


d(m)d(n)d(n  —  m) 


0,  \m\  >  L 


(4.14) 


(4.15) 


The  frequency  and  time  domain  responses  of  the  two  dimensional  window  are  shown 
in  Figures  4. 1-4.2.  Notice  that  the  optimum  window  has  a  wide  main  lobe  and  low 
sidelobes. 


4.3.2  Direct  Method 

Compute  the  impulse  response  T(r): 


i  N- 1 

Y(r)  =  —  X(fc) exp (—  j4nkr/cN) 
™  k= o 


(4.16) 


A  Hanning  window  d(m)  may  be  applied  to  the  frequency  samples 

,,  .  1,  2n(m  —  (N  —  l)/2) 

d(m)  =  -(1  +  cos( — - - jX 


(4.17) 


then  the  bispectrum  can  be  computed  over  the  triangle  0  <  rx  <  r2,  rj  +  r2  <  fs/ 2 
where  f,  is  the  sampling  frequency. 
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B{ru  r 2)  =  Ky{rx)y{r2)y\ri  4-  r2)  (4.18) 

where  I\  is  constant  and  depends  on  the  spacing  between  frequency  samples  (K  = 
( j-)2),  fs  is  the  sampling  frequency  of  the  data. 

For  sufficiently  large  number  of  data  points,  both  methods  provide  unbiased  es¬ 
timates.  In  general  classical  methods  have  high  variance  and  the  estimates  obtained 
are  not  usually  smooth.  The  advantage  of  using  classical  techniques  is  the  fact  that 
fast  computations  algorithms  such  as  the  FFT  are  available. 

4.4  Blade-Sphere  Example 

Consider  a  blade  of  width  rx  —  r0  and  a  sphere  at  a  distance  r2  —  rx  from  the 
trailing  edge  of  the  blade.  The  zero  reference  is  set  at  a  distance  r0  from  the  leading 
edge  of  the  blade  (see  Figure  4.3).  It  can  be  shown  that  if  one  considers  four  ray 
paths  the  backscattered  signal  is; 

St(fj)  =  Ai  +  A2exp(-j27rr1^-)  (4.19) 

c 

f. 

-t-A3exp(-j27rr2-i) 

c 

+A4  exp  (~j2x(rl  +  r2  -  r0)— ) 

c 

assuming  that  the  amplitudes  Ax,  A2,  A3,  and  A4  are  frequency  independent.  The 
corresponding  amplitude  and  phase  of  this  structure  are  shown  in  Figure  4.4.  The 
impulse  response  of  the  blade  and  sphere  target  as  computed  using  the  inverse 
Fourier  transform  (with  a  Hanning  window  over  a  frequency  band  from  1.5  GHz 
to  12  GHz)  is  shown  in  Figure  4.5.  Note  the  one  to  one  relationship  for  the  first 
three  responses  with  the  geometry  of  the  target.  The  fourth  response,  however,  is  an 
interaction  term  and  there  is  no  single  discrete  target  component  which  corresponds 
to  it. 

Figure  4.6  shows  the  bispectral  response  of  the  backscatter  from  the  blade-sphere 
combination  without  using  a  window.  The  bispectral  response  of  the  windowed  third 
moment  sequence  is  shown  in  Figure  4.7.  In  these  figures,  r0  =  0  (zero  reference  at 
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leading  edge  of  the  blade).  Note  the  response  at  (ri  =  1.6meters,  r2  =  2.4meters) 
which  is  due  to  interaction  between  the  trailing  edge  of  the  blade  and  the  sphere 
and  is  represented  by  the  term  exp(— j2-K(rx  +  r2  —  r0)^-).  We  can  also  notice  the 
twelve  symmetry  regions  of  the  bispectrum. 

If  the  value  of  r0  is  increased  to  1  meter,  the  bispectral  response  shown  in  Figure 
4.8  results.  Note  that  all  of  the  response  terms  have  moved  along  the  associated 
diagonals,  but  the  particular  scattering  and  interaction  terms  are  still  observable. 
A  three  dimensional  display  of  the  bispectral  response  is  shown  in  Figure  4.9. 

4.5  Bispectral  Response  of  Noisy  Data 

If  zero  mean  white  Gaussian  noise  is  added  to  the  raw  data,  the  resulting  impulse 
response  and  bispectral  response  are  shown  in  Figures  4.10  and  4.11  respectively. 
Note  that  the  detectability  of  the  target  scattering  terms  is  greatly  reduced  in  the 
impulse  response,  and  noise  components  in  the  time  domain  can  be  misinterpreted 
as  discrete  scatterers.  This  problem  is  less  effective  in  the  bispectrum  case  where 
discrete  scatterers  and  interaction  terms  can  still  be  identified  under  the  same  noise 
conditions  used  in  the  impulse  response.  This  is  due  to  the  fact  that  in  the  bispec¬ 
trum  case  the  noise  components  at  different  frequencies  are  uncorrelated. 

4.5.1  Mean  Square  Error 

Consider  the  measured  observation  y(f)  at  frequency  /.  Assume  that  the  noisy 
data  is  y(f )  =  x(f)  +  nj  where  nj  is  a  Gaussian  distributed  random  variable.  Also, 
assume  that  the  data  are  wide  sense  stationary.  The  frequency  bandwidth  is  fh  —  fi 
where  //  >  0  Hz. 

The  third  order  autocorrelation  function  of  the  measured  noisy  data  is; 

Ryifu  fi)  =  E{y(f)y(f  +  h)y{f  +  /2)}  (4.20) 

then 

Ry'fuh)  =  E{x(f)x(f  +  /!)*(/  +  f2) 

+  X(/)n/  +  /in/+/2  +  X(/  +  f\)ninf+h 
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+  x(f  +  f2)n}nJ+h  +  x{f)x(f  +  fi)nf+h 
+  x(/)x(/  +  /2)n/+/l  +  x(f  +  fi)x(f  +  f2)nf 
+  nlnf+hnl+h 


(4.21) 


therefore 

Ry(fl,  h)  =  RX(fu  h)  +  HRnifx)  +  Rnifi)  +  Rn(fl  ~  /2)}  +  0Rn(fl,  /j)(4.22) 

However,  for  a  Gaussian  noise  0  =  0  then 

/2)  =  J?,(/i,  /2)  +  <r2x{*(/i)  +  *(/2)  +  6(f 1  -  /2)}  (4.23) 

where  <r2  is  the  noise  power. 

The  second  order  autocorrelation  is  given  as 

Ry(F)  =  E{y(f)y(f  +  F)}  (4.24) 

=  E{(x(f)  +  nf)(x(f  +  F)  +  nf+F)} 

=  RX{F)  +  Rn(F) 

=  RX(F)  +  a26(F) 

Consider  the  following  as  an  error  measure  for  estimating  the  impulse  response  for 
noisy  targets; 

Ej  =  E{Jt  | K(t)  -  hy(t)\2dt}  (4.25) 

where  hy(t)  and  hT(t)  represent  the  impulse  response  for  noisy  and  noiseless  data 
respectively  and  E{.}  stands  for  the  expected  value.  Then  Ej  can  be  written  as  ; 

E,  =  £{  jf  \H,U)  -  Hu(J)\2df }  (4.26) 

using  Parseval’s  theorem,  then  E\  is  equal  to; 

E,  =  E{  /'*  \n(f)\2df)  (4.27) 

J  h 

Ei  =  2  a2(fk-ft) 
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If  the  spectral  density  is  computed  then  the  error  measure  is; 

Es  =  Jt\Sx(t)  -  Sy(t)\2dt 

=  I  \Rx(f)  -  Ry(f)M 

=  fys(f)  m 

=  0 


Notice  that  in  this  case  the  mean-square  error  is  not  dependent  on  the  bandwidth. 
Consider  the  following  as  an  error  measure  for  bispectrum  of  noisy  signals; 

Eg  =  v/i  X  "  Ey(ti,t2)\2dtidt2  (4.28) 

where  Bx(t\,  <2),  By(ti,  1 2)  axe  the  bispectrum  for  ideal  noiseless  data  and  noisy  data 
respectively.  Then  using  Paxseval’s  theorem  in  two  dimensions  Eb  can  be  written 
as; 

Eb  =  (4.29) 

Then  using  equation  (4) 

Eb  =  jf  jh  |<(/,)  +  <(/,)  +  «(/.  -  ft) MMh  (4.30) 

=  2<t!xv'(2  +  V2)(/»  -  /1)  (4.31) 


Therefore  : 


Eb  =  Ejx 


2  +  V2 

M  h-fi 


(4.32) 


Therefore  the  mean  square  error  is  directly  dependent  on  the  frequency  band¬ 
width  in  the  impulse  response  case  as  well  as  the  bispectral  response.  However,  in 
the  bispectrum  case  a  new  factor  is  included  that  depends  on  the  signal  itself. 
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4.6  Bispectral  Response  For  Dispersive  Scatterers 


The  effect  of  dispersion  in  the  data  is  shown  in  Figures  4.12  and  4.13.  Disper¬ 
sion  is  the  result  of  scattering  from  a  subcomponent  of  the  radar  target  where  the 
scatterer  behavior  varies  as  a  function  of  frequency.  In  this  case,  the  behavior  is 
modeled  as  follows. 

Ai(fj)  =  Ai(l  +  a  exp  (  —  kfj)  cos(/3fj))  (4.33) 


where 

Ai{fj )  =  is  the  frequency  dependent  amplitude  of  the  backscatter  coefficient 
from  the  ith  scatterer.  and 
a,  A*,  and  0  are  constants. 

In  this  study,  a  =  0.6  and  k  =  0.01,  and  0  =  207T. 

These  figures  show  that  the  effect  of  dispersion  is  not  critical  to  the  bispectrum 
response. 

4.6.1  Single  Edge  Dispersion 

The  above  form  of  dispersion  effects  the  amplitude  of  the  response  only.  In 
general  dispersion  effects  the  amplitude  and  the  phase  of  the  radar  back-scatter. 
This  type  of  dispersion  results  in  energy  loss  along  the  physical  target  thus  effectively 
moving  the  scattering  centers  in  the  down  range  profile. 

In  this  section  we  study  the  case  where  only  one  edge  is  dispersive  in  amplitude 
and  phase.  This  is  rather  a  realistic  model  and  matches  to  a  certain  extent  scattering 
from  a  rocket  where  the  front  end  scatterer  (the  nose)  is  highly  dispersive.  In  this 
case  the  amplitude  of  the  returned  signal  can  be  modeled  as: 

Mfj)  =  Ai{(l  +  a  exp  (-£/,)  cos( /?/,)}  exp (jp/f)  (4.34) 

where  p  depends  on  the  scatterer.  This  model  indicates  that  scattering  takes  place 
near  the  point  of  the  nose  for  high  frequency  components  and  further  away  from 
that  point  at  low  frequencies.  Figures  4.14-4.15  show  the  impulse  response  and  the 
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bispectral  response  of  this  type  of  scattering.  As  expected  the  scatterers  location 
will  be  displaced  and  the  response  will  be  smeared.  The  bispectral  response  shows  a 
similar  type  of  behavior  except  that  the  energy  loss  phenomenon  and  the  frequency 
dependent  scattering  is  recognizable.  Notice  that  the  interaction  term  response 
inuicaies  uhat  scat  wring  place  somewhere  between  the  nose  and  the  second 

scatterer. 

4.7  Bispectral  Responses  for  Canonical  Aircraft 

Experimental  data  for  scattering  from  a  6  Inches  canonical  aircraft  for  frequen¬ 
cies  between  1  —  18  GHz  were  used  in  determining  the  bispectral  responses  from 
various  components  of  the  aircraft.  Back-scattering  from  wing  only,  fuselage,  fuse¬ 
lage  and  stabilizer,  and  complete  aircraft  parts  were  measured  at  the  Ohio  State 
ElectroScience  laboratory  [13]. 

Figures  4.16-4.19  show  the  bispectral  response  for  such  models  corresponding  to 
different  parts  of  the  aircraft.  In  the  case  where  the  scattering  from  the  fuselage 
only  is  measured  no  considerable  interaction  term  is  observed.  The  addition  of  the 
stabilizer  introduces  some  interactions.  In  the  case  where  only  the  wing  is  present,  it 
is  shown  that  many  high  order  scattering  mechanisms  take  place  along  the  wing  area. 
The  bispectral  response  for  the  complete  aircraft  is  dominated  by  the  interactions 
that  take  place  in  the  wing  area.  These  figures  will  help  us  understand  bispectral 
responses  from  real  targets  where  many  complicated  scattering  mechanisms  are 
involved. 

4.8  Bispectral  Responses  for  Commercial  Aircraft 

The  compact  radar  cross  section  measurement  range  at  the  Ohio  State  Univer¬ 
sity  ElectroScience  Laboratory[13]  was  used  to  make  measurements  of  the  radar 
scattering  of  a  set  of  5  scale  models  of  commercial  transport  aircraft  (see  Figure 
4.20).  The  measurements  were  calibrated  so  that  absolute  values  of  amplitude  (in 
square  centimeters)  and  phase  (in  degrees  from  an  absolute  reference  plane)  were 
available  over  the  band  of  frequencies  from  1.5  to  12.0  GHz.  This  corresponds  to 
frequencies  in  the  HF  band  for  full  scale  aircraft  and  means  that  the  targets  are  in 
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the  resonance  region  of  the  radar  (ie:  The  wavelength  is  on  the  order  of  the  target 
size)[13]. 

Using  the  techniques  described  above  (including  the  use  of  the  optimum  win¬ 
dow  function),  the  bispectrum  of  each  radar  target  signature  was  computed.  The 
results  are  n'uowu  in  Figures  4.21  to  4.25.  Note  that  the  different  radar  targets  yield 
distinctly  different  bispectral  responses.  It  can  be  shown  that  this  implies  that  such 
responses  are  likely  candidates  as  features  for  radar  target  identification. 

Since  each  radar  target  is  being  measured  in  the  resonance  region  (as  mentioned 
above),  it  can  be  shown  that  the  target  signature  is  relatively  insensitive  to  changes 
in  aspect  angle,  (especially  the  impulse  response  signature.)  The  result  of  taking 
the  bispectral  response  for  the  Boeing  707  as  the  aspect  angle  varies  from  10  to  30 
degrees  is  shown  in  Figures  4.26-4.28.  It  can  be  seen  that  the  bispectral  response  is 
rather  stable  between  each  10  degree  angle  change.  Figure  4.29  shows  the  bispectral 
response  of  a  747  at  0  degrees  azimuth  and  -5  dB  signal  to  noise  ratio  (SNR).  In 
the  next  section  bispectral  images  (as  such)  are  used  in  classifying  different  targets. 

4.9  Classification  Methods 

In  this  section  we  describe  how  to  classify  radar  targets  based  on  their  bispectral 
response  using  two  dimensional  cross-correlation.  The  concept  of  the  “catalog”  of 
data  will  be  used  with  a  bispectral  feature  space.  In  this  technique  a  catalogue  of 
noise  free  bispectral  responses  are  used  to  determine  the  identity  of  a  target  from 
which  a  noisy  bispectrum  is  available.  It  is  known  (in  this  study)  that  the  target 
corresponds  to  one  of  the  entries  in  the  catalogue.  First,  the  goal  to  find  a  catalogue 
bispectral  response  that  most  closely  matches  the  bispectral  response  BJ(rl,r2)  of 
an  unknown  target.  That  is  to  minimize 


min {JT  J  r2)  -  BJ(rl,r2))2dr1dr2  = 

(4.35) 

min  {J  J  Bf(rur2)drldr2  +  J  J  B]{r1,r2)dr1dr2 

(4.36) 

~  J  J  Bt{rur2)Bj{rl,r2)drldr2} 

(4.37) 
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where  Bi(ri,r2)  and  Bj(r i,r2)  axe  the  bispectral  responses  of  catalogue  target  i  and 
unknown  target  j.  This  requires  maximizing 


I  [  Bi(rl,r2)Bj(r1,r2)drldr2 

Jr>  J  to 


(4.38) 


and  *hic'  fa p  he  done  by  shifting  Pj(r'i,r2)  over  all  the  domain  of  B.fr j  r2)  This  is 
equivalent  to  maximizing 

In  It,  Bi(rur2)Bj(ri  +  Rur2  +  R2)drxdr2 


Pij(R  1)  -^2)  — 


[/,,  fr,  |Bj(n  +  Rur,  +  Rt)Vdndr,\'  [/„,  /„  |Bi(r„rJ)|yr1*J 


-(4.39) 


where  pij(R\,  R2)  is  the  normalized  cross  correlation  coefficient  of  catalogue  target 
bispectral  response  i  =  1, 2, . . . ,  5  and  test  target  response  j.  Notice  that  the  factor 
in  the  denominator  depends  on  R\,  R2  and  this  is  in  order  to  compute  a  normalized 
crosscorrelation. 


4.9.1  Indirect  method 

Using  Fourier  transform  identities  and  Parsevals  theorem  the  cross  correlation 
can  be  written  as 


PlJ(RuR,)  = 


[//.  4  mh,h)Vdhdhf  [4  4  \Ri(fiJ2)M,df2}i 


(4.40) 


where  I2FT  stands  for  the  inverse  two  dimensional  Fourier  transform  and  Rj(fi,  f2) 
is  the  conjugate  of  Rj(fi,  /2).  For  the  discrete  frequency  case  the  cross  correlation 
is 

I2FFT{RiUuh)R'(h,f2)} 


Ptji.R\i  R2) 


[ea  e*  mh,f2)V\’  [E/,  eab>(/../>)ij]} 


(4.41) 


where  I2FFt  is  the  inverse  two  dimensional  fast  Fourier  transform.  The  test  target 
is  classified  to  a  catalogue  class  c  if 


pcj(RuR2)  =  max{max  pij(Ru  R2)}  i  =  1,...,5. 

•  Hi  ,«2 


(4.42) 
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4.9.2  Direct  method 


If  the  direct  method  for  estimating  the  birange  profile  for  radar  targets  was  used 
then  the  classification  would  follow  as 

•  Compute  the  birange  estimate  of  both  catalogue  target  i  and  test  target  j  as 


x?,(rbr2)  =  /3y,(r1)yl(r2)yi(r1  +  r2) 


(4.43) 


over  the  triangle  A  :  jq  <  r2,  rq  +  r2  <  cn,  and  r2  >  0.  where  /3  is  a  constant 
and  y(r)  is  the  impulse  response, 


Vi{r )  =  -rTY.  Xi(n)exp(j2nnr/N) 

n= 0 


(4.44) 


and  X(n)  is  the  measured  frequency  data  for  target  i. 

Compute  the  crosscorrelation  p,j(i?i,722)  as 

(rt  O  X _ Er,  £r2  r2)BJ(r1  +  Ru  r2  +  R2) _ 

"  15  2  EnEr2|Bi(r1,r2)P]MEnErJBj(rl  +  i?1,r2  +  H2)Pl^  ^ 

where  Rv,  R2  C  A  . 


•  The  decision  rule  is; 


pCJ(RuR2)  =  max{max  -R2)}  i  =  1,...,5. 

t  ri\  ,I\2 


(4.46) 


4.9.3  Comparison 


The  first  method  does  not  require  bispectrum  estimation  but  still  requires  esti¬ 
mates  of  third  order  autocorrelation,  while  the  second  method  requires  computation 
of  two  dimensional  crosscorrelation  but  no  third  order  autocorrelation.  Also,  the 
first  method  involves  a  two  dimensional  window,  while  second  method  involves  a 
bispectral  estimate  of  Hanning  window.  Another  important  difference  is  that  the 
first  method  requires  computation  of  third  order  autocorrelation  over  all  birange 
domain  while  second  method  makes  use  of  the  triangular  symmetry.  Moreover, 
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the  first  method  requires  third  order  correlation  over  all  iZj,  -R2,  two  dimensional 
Fourier  transform.  The,  second  method  requires  one  dimensional  Fourier  transform, 
birange  estimate  over  a  triangle  A,  and  crosscorrelation  estimate  of  biranges  over 
A. 

In  order  to  compute  a  normalized  cross  correlation  coefficient  we  need  to  com¬ 
pute  the  normalizing  factor  Vi?i,i?2  6  A.  This  is  a  major  disadvantage  for  the 
second  method.  Also,  the  decision  criteria  for  the  first  method  is  over  all 
while  the  other  method  is  over  one  triangle(l/8)  of  domain. 

4.9.4  Classification  Algorithm 

It  has  been  shown  that  both  classical  methods  in  bispectrum  estimation  require 
many  computations  if  used  for  radar  target  identification.  An  algorithm  that  com¬ 
bines  both  approaches  can  be  implemented  as  follows; 

1)  Compute  the  Impulse  response  of  the  test  target  j  using  discrete  Fourier  trans¬ 
form  which  can  be  implemented  as  Fast  Fourier  transform. 

i  N~} 

Zj(r)  =  —  Xj(f)  exo  j4nfr/c - FFT  (4.47) 

j= 0 

2)  Compute  the  third  order  autocorrelation  estimate 

RjUufo)  =  -^X^Mri)*;(r2)Mri  +  r2)expj47r(/1r1  +  /2r2)/c - 2FFT^ .  48) 

3)  Compute  the  cross  correlation  coefficient  as 

[E;,  mfx,h)V\k  [£,,  ^:hR,(h.h)\2f 

4)  The  test  target  is  classified  to  a  catalogue  class  c  if 

Pcj{R\,  Ri)  =  max{max  /?2)}  i  =  l,...,5.  (4.50) 

*  R\  ,/?2 

In  the  above  algorithm  one  one-dimensional  FFT  and  two  two-dimensional  FFt 
are  needed  in  order  to  compute  the  normalized  cross  correlation  coefficient.  Thus, 
the  above  algorithm  is  most  suitable  for  classification  of  large  arrays  of  data  in  the 


2FFT  (4.49) 
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frequency  domain  where  the  computation  time  saved  due  to  the  parallel  computa¬ 
tion  feature  of  the  above  algorithm  becomes  considerably  large. 

4.10  Radar  Target  Identification  Example 

A  Monte  Cailo  simulation  based  on  the  Indirect  Method  was  performed  to  eval¬ 
uate  the  performance  of  the  concepts  discussed  in  this  chapter  where  the  bispectral 
response  was  used  as  a  feature  for  characterizing  radar  targets.  The  data  base  used 
in  this  example  has  been  used  before  in  radar  target  identification  studies  [14,15,16]. 
The  data  base  consists  mainly  of  experimental  measurements  in  the  frequency  band 
from  1  —  12  GHz  of  scale  models  of  commercial  transport  aircraft.  A  photo  with 
corresponding  names  of  these  aircraft  are  shown  in  Figure  4.20.  Moreover,  the  scale 
data  corresponds  to  measurements  of  the  radar  cross  section  (RCS)  of  full  scale 
aircraft  in  the  HF /VHF  frequency  band. 

The  classification  algorithm  is  shown  in  Figure  4.30.  In  this  algorithm  Gaussian 
distributed  noise  was  added  to  the  raw  data  (complex  radar  cross  section  mea¬ 
surements  RCS  versus  frequency).  An  estimate  of  the  bispectral  response  of  the 
resulting  noisy  data  is  then  computed  using  the  indirect  method.  A  two  dimen¬ 
sional  cross-correlation  was  performed  between  the  noisy  bispectral  response  and 
the  catalogue  ideal  bispectral  response.  A  decision  is  made  in  favor  of  the  catalogue 
with  the  largest  cross  correlation.  The  same  procedure  is  repeated  for  different 
test  targets.  Figure  4.31  shows  the  actual  algorithm  used  to  evaluate  the  classi¬ 
fication  performance  technique  shown  in  Figure  4.30.  In  this  algorithm  the  third 
moment  sequence  of  both  noise  free  catalog  data  noisy  data  (additive  noise)  are 
computed.  The  normalized  inverse  two  dimensional  fast  Fourier  transform  of  the 
product  of  the  two  third  order  moment  sequences  is  then  computed.  Decision  is 
made  in  favor  of  the  class  with  the  largest  cross  correlation  with  the  unknown  noisy 
target.  Decision  statistics  for  each  target  are  computed  at  each  noise  level  and  total 
statistics  of  classification  error  are  computed  for  all  targets  at  a  certain  noise  level. 
One  hundred  experiments  were  performed  for  each  test  target  (total  of  500  experi¬ 
ments).  The  entire  test  is  repeated  for  different  noise  levels.  Finally  the  percentage 
of  misclassification  is  plotted  versus  signal  to  noise  ratio. 


44 


Figure  4.32  shows  the  classification  performance  for  five  commercial  aircraft 
where  complete  azimuth  information  is  assumed  to  be  known.  This  performance 
curve  shows  that  bispectral  estimates  can  be  used  effectively  in  radar  target  identi¬ 
fication.  Notice  that  the  percentage  of  misclassification  is  zero  for  a  signal  to  noise 
higher  than  0  dB. 

The  above  example  shows  that  reliable  target  identification  can  be  achieved 
using  the  bispectral  response  as  a  feature  source.  The  classification  performance 
can  be  further  improved  by  developing  a  feature  selection  algorithm  that  can  be 
applied  to  the  bispectral  responses  of  both  catalogue  and  unknown  test  targets. 
Furthermore,  the  classification  approach  used  in  this  example  is  a  non-parametric 
technique  that  requires  no  a  priori  information  about  the  parametric  distribution 
of  the  measured  radar  data. 

4.11  Conclusions 

In  this  chapter  the  application  of  bispectrum  estimation  to  radar  signature  anal¬ 
ysis  is  investigated.  It  has  been  demonstrated  that  a  bispectral  response  of  radar 
backseat  ter  signals  has  the  advantage  of  identifying  multiple  interactions  for  com¬ 
plex  radar  targets.  This  leads  to  the  conclusion  that  higher  order  interactions  of 
radar  backscatter  may  be  identified  using  higher  order  spectral  analysis.  Such  high 
order  interactions  that  are  relatively  smaller  in  the  impulse  response  because  of  their 
amplitudes  may  be  recovered  using  high  orders  of  spectral  analysis.  Therefore,  the 
bispectrum  and  other  high  order  spectral  analysis  methods  may  be  considered  as 
potential  radar  feature  extraction  techniques  that  may  be  effectively  used  jointly 
with  the  impulse  response  in  radar  target  identification.  It  is  shown  that  bispectral 
responses  of  radar  backscatter  signals  are  relatively  robust  under  conditions  of  high 
noise  and  frequency  dispersion.  Finally,  it  has  been  shown  that  radar  targets  can  be 
successfully  identified  at  zero  dB  signal  to  noise  ratio  using  a  simple  classification 
algorithm  such  as  cross  correlation. 
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Figure  4.1:  Optimum  window  (minimum  bias  supremum)  in  frequency  domain 
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Figure  4.3:  Internal  bounce  diagram  for  blade-sphere  interaction. 
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Figure  4.4:  Frequency  response  for  blade-sphere  target. 
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Figure  4.6:  Bispectrum  for  blade-sphere  without  window. 
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Figure  4.7:  Bispectrum  for  blade-sphere  r\  =  1.6  m,  r2  =  2.4  m.  with  window 
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Figure  4.8:  Bispectrum  for  blade-sphere  with  r0  = 
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Figure  4.10:  Impulse  response  for  noisy  backscatter  from  blade-sphere  (SNR=  -30 
dB) 
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Figure  4.11:  Bispectrum  for  blade-sphere  with  noisy  data  (SNR=  -30  dB) 
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Figure  4.12:  Impulse  response  of  frequency  dispersive  backscatter  for  blade-sphere. 
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Figure  4.13:  Bispectrum  of  frequency  dispersive  backscatter  for  blade-sphere 
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Figure  4.15:  Bispectrum  of  frequency  dispersive  backscatter  for  blade-sphere  with 
single  edge  dispersion 
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Figure  4.16:  Bispectral  response  for  fuselage  at  0  degrees  azimuth,  0  degrees  eleva¬ 
tion. 
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Figure  4.18:  Bispectral  response  for  wing  at  0  degrees  azimuth,  0  degrees  elevation. 
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Figure  4.19:  Bispectral  response  for  complete  canonical  aircraft  at  0  deg.  azimu 
0  deg.  elevation. 
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Figure  4.20:  Scaled  drawings  for  aircraft  used  in  experimental  study 
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Figure  4.21:  Bispectral  response  for  747  at  0  degrees  azimuth,  0  degrees  elevation 
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Figure  4.23:  Bispectral  response  for  concord  at  0  degrees  azimuth,  0  degrees  eleva¬ 
tion. 
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Figure  4.24:  Bispectral  response  for  DC10  at  0  degree  azimuth,  0  degree  elevation. 
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Figure  4.25:  Bispectral  response  for  707  at  0  degrees  azimuth,  0  degrees  elevation. 
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26:  Bispectral  response  for  707  at  10  degrees  azimuth,  0  degrees  elevation 
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Figure  4.27:  Bispectral  response  for  707  at  20  degree  azimuth,  0  degree  elevation. 
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Figure  4.28:  Bispectral  response  for  707  at  30  degree  azimuth,  0  degree  elevation 
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Figure  4.29:  Bispectral  response  for  707  at  0  degrees  azimuth,  0  degrees  cleva 
and  SN  R  =  —5 dB. 
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Figure  4.31:  Diagram  of  bisprctral  radar  target  classification  algorithm  as  iui|>l< 
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Figure  4.32:  Classification  performance  assuming  complete  azimuth  information. 
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5.  Investigations  of  Structural  Pattern  Recognition  Methods  for 

Radar  Target  Identification 


5.1  Introduction 

One  of  the  more  troublesome  aspects  of  the  level-crossing  based  preprocessing 
methods  of  [17]  is  the  determination  of  conditional  string  probabilities  by  the  his¬ 
togram  method.  Convergence  of  the  histograms  as  well  as  the  issue  of  unclassifiable 
strings  unnecessarily  complicates  the  investigation  of  inference  of  more  sophisti¬ 
cated  language  theoretic  classifiers  and  the  computation  of  the  resulting  classifiers 
performance.  To  circumvent  these  difficulties,  analytical  methods  for  determination 
of  unsqueezed  octant  crossing  and  single  level  crossing  have  been  developed. 

Currently,  the  algorithms  are  being  used  to  determine  misclassification  prob¬ 
abilities  by  Monte-Carlo  simulation.  This  method  potentially  allows  for  shorter 
execution  time  than  the  histogram  methods  of  likelihood  function  estimation  used 
previously,  since  conditional  string  probabilities  need  only  be  calculated  for  those 
strings  which  occur  during  *he  Monte-Carlo  test.  Furthermore,  the  problem  of 
un- classifiable  strings,  as  well  as  the  problems  related  to  the  matched-noise  phe¬ 
nomenon  are  alleviated.  The  algorithms  which  are  employed,  their  current  and 
possible  future  uses  are  described  in  the  following  section. 

While  extremely  compact,  the  non-deterministic  form  of  the  a  syntactic  classifier 
derived  in  [17]  is  not  particularly  revealing  in  terms  of  analyzing  decision  rules 
induced  by  such  a  classifier.  In  order  to  gain  insight  into  the  resulting  decision  rules, 
programs  have  been  developed  which  convert  these  classifiers  to  a  deterministic 
form.  The  deterministic  form  resembles  a  tree  classifip’’  4s  such,  it  reveals  the 
relative  discriminating  power  of  particular  symbol  ^strings.  This  notion  has 

yet  to  be  formalized. 

The  deterministic  form  of  the  classifier  may  be  used  as  an  analysis  tool  to  help 
understand  how  the  grammatical  inference  methods  are  able  to  extend  decision 
rules  generated  by  likelihood  function  estimation.  Recall  that  when  likelihood  func¬ 
tions  estimated  by  the  histogram  method  are  used  for  decision  rule  generation  (i.o. 
decisions  do  not  exist  for  every  possible  symbolic  representation)  the  grammati- 
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cal  inference  method  allows  the  decision  regions  to  expand  into  areas  of  unknown 
classification  to  reduce  computational  complexity  (see  [17]).  Inspection  of  the 
deterministic  form  of  the  classifier  reveals  the  extent  to  which  overlapping  of  the 
decision  regions  induced  by  the  grammatical  inference  method  takes  place.  When 
coupled  with  the  analytical  method  of  likelihood  function  estimation,  this  form  of 
the  classifier  will  be  able  to  quantify  the  gains  realized  by  the  expanded  classifier. 
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Figure  5.1:  Conditional  Probability  Overlaps 

5.2  Analytical  Methods  for  Likelihood  Function  Determination 

The  basic  algorithm  is  designed  to  calculate  the  conditional  string  probability 
for  each  element  of  a  catalog  to  within  a  given  tolerance.  The  string  probability 
calculation  procedures  use  this  specified  tolerance  to  determine  allowable  error  for 
each  calculation  in  the  process.  These  individual  errors  are  set  so  that  the  total 
string  probability  estimate  is  within  the  specified  tolerance.  The  underlying  pur¬ 
pose  for  calculation  of  the  conditional  string  probabilities  is  to  make  a  decision  as 
to  which  catalog  element  the  string  was  generated  under.  The  maximum  likelihood 
estimate  of  the  generating  catalog  element  is  determined  by  choosing  the  maximum 
of  the  conditional  probabilities.  A  decision  can  only  be  made  if  the  band  about  the 
iriaximu  n  conditional  probability  (the  width  determined  by  the  given  conditional 
probability  tolerance)  does  not  overlap  any  other  bands  (see  Figure  5.1).  Since  little 
is  known  about  a  given  conditional  string  probability  before  estimation,  the  total 
string  probability  tolerance  is  first  set  to  some  nominal  value  and  the  conditional 
string  probabilities  are  determined.  If  the  situation  exists,  where  the  maximum 
conditional  probability  overlaps  some  other  conditional  range  of  probabilities  then 
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the  conditional  probabilities  are  re-calculated  with  an  appropriately  decreased  tol¬ 
erance.  In  this  way,  the  conditional  string  probability  estimation  algorithms  may 
be  used  to  make  a  Bayesian  decision  about  catalog  elements. 

Clearly,  the  string  probability  calculation  algorithms  may  be  used  for  Monte- 
Carlo  test  since  it  allows  for  a  hard  decision  of  any  given  string.  However,  these 
algorithms  will  allow  for  grammatical  inference  with  supervision,  since  the  string 
probability  routines  may  act  as  a  teacher.  Furthermore,  since  the  conditional  proba¬ 
bilities  are  available  on  demand,  investigation  of  information  content  of  the  symbolic 
representation  is  allowed  for  as  well  as  investigation  of  the  precise  influence  of  noise 
on  the  likelihood  functions. 

5.2.1  Conditional  String  Probability  Estimation 

Each  of  the  probability  calculation  algorithms  involve  further  layers  of  error- 
control  iteration  beyond  the  decision  process  described  in  the  section  5.2.  This 
sub-section  describes  the  first  two  of  these  layers.  They  are  common  to  both  the 
octant  crossing  and  single  level  crossing  schemes  and,  with  modification,  would 
apply  to  any  symbol  assignment  scheme  with  adaptive  thresholds. 

The  first  layer  of  error  control  iteration,  accounts  for  variations  in  the  magnitude 
threshold  by  approximating  it  with  a  normal  random  variable.  The  thresholds  of 
the  double  level  crossing  are  dependent  upon  the  maximum  and  the  minimum  of 
observed  measurements  and,  as  such,  involve  order  statistics  of  independent  Rician 
random  variables. 

The  second  layer  of  error  control  iteration  involves  the  relation  between  individ¬ 
ual  symbol  error  and  total  string  error.  The  technique  developed  in  this  section  is 
valid  for  any  symbol  assignment  scheme  in  which  each  symbol  corresponds  to  an 
independent  event.  The  individual  symbol  probability  estimation  procedures  are 
described  in  the  following  section. 

Magnitude  Threshold  as  a  Random  Variable 

By  the  definition  of  the  single  level  crossing  and  octant  crossing  schemes,  the 
magnitude  threshold  used  in  these  schemes  is  the  sum  of  the  individual  magnitudes 
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Equation  (5.1).  Since  the  threshold  used  for  string  generation  is  a  function  of  the 
observation,  it  too  is  a  random  variable.  Thus  conditional  string  probabilities  are 
calculated  under  the  assumption  of  a  fixed  threshold  and  integrated  with  respect 
to  the  threshold  density  to  determine  the  probability  of  a  given  string  conditioned 
only  on  catalog  element. 

1  v 

r=vZK«‘)l  (5-D 

*V  !=t 

In  equation  (5.1),  T  is  the  magnitude  threshold  N  represents  the  number  of 
measurements  taken  and  a,  is  the  ith  noisy  measurement. 

P(S  |  Q)  =  /X  P(S  |  T  =  f, C,)P(T  =  t  |  Ci)dt  (5.2) 

The  individual  magnitudes  of  the  observation  vector  are  Rician  distributed  ran¬ 
dom  variables,  by  assumption  of  independent  Gaussian  noise.  Equation  (5.1)  in¬ 
dicates  that  the  magnitude  threshold  is  a  weighted  sum  of  these  Rician  random 
variables.  Thus,  by  the  central  limit  theorem,  it  is  reasonable  to  assume  that  the 
magnitude  threshold  for  both  schemes  is,  conditionally,  a  Gaussian  distributed  ran¬ 
dom  variable.  Furthermore,  the  central  limit  theorem  indicates  the  statistics  of  the 
magnitude  threshold. 


P(T  =  t  |  Cj)  = 


(5.3) 


And  so  Equation  (5.2)  becomes: 
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P(S  I  C\)  = 


y/2tt<7? 


/OO 

P(S 
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t,  Ct)e 


■M' 


dt 


(5.4) 


In  equations  (5.3)  and  (5.4)  jit  i  -  the  threshold  mean  and  a?  is  the  threshold  variance. 
They  are  calculated  from  components  of  the  measurement  vector  mean  ,E  [c,2],  and 
second  moment,  E[zt]2. 


dt 
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,v 


(5.5) 
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I  —  1 


(5.6) 


The  components  of  the  measurement  vectors’  mean  and  second  moment  are,  in 

turn  functions  the  magnitude  of  components  of  the  noiseless  measurement  vectors, 

2 

fj.,t .  and  the  per-channel  noise  variance,  a2  (the  quantity  g,  is  the  ratio  jf?)- 


E  k] 


hz,  +  2(7 


2 


(5.7) 

(5.8) 


The  Gaussian  Quadrature  method  is  used  to  approximate  the  ir‘?gral  of  Equa¬ 
tion  (5.2)  with  a  sum  (see  Reference  [18]). 

m 

P(S  I  Ci)  «  £ff*P(S  I  T  =  *k,Ci)P{T  =  xk  i  Ci)  (5.9) 

fc=  i 

This  method  chooses  the  abscissas  (xk)  and  the  weights  ( Hk )  in  such  a  way  as  to 
minimize  the  error  between  the  integral  and  the  sum;  exact  answers  are  obtained  for 
polynomials  of  degree  m  —  1  or  less.  Gaussian  Quadrature  with  Hermite  weighting 
is  used  to  approximate  integrals  with  weighting  functions  of  the  form  iv(t)  —  e~  '2. 
Thus  a  change  of  variable  is  made  to  put  Equation  (5.4)  into  such  a  form. 


x  =  h(t) 


t  ~  fh 

°t 


(5.10) 


t  —  h  x(x)  =  atx  -f  fit 


(5.11) 


And  so  Equation  (5.4)  becomes: 

P(S  I  Ci)  =  ~  r  P(S  I  T  =  h-\x),Ct)e-x2dx  (5.1 2) 

v2tt  9-oo 

Gaussieui  Quadrature  with  Hermite  weighting  may  now  be  used  to  approximate 
the  integral  of  Equation  (5.12): 

m 

P(S  I  C\)  «  £  HkP(S  I  T  -  h-'(xk),  Ci)  (5.13) 

k=  1 
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The  error  associated  with  an  estimate  using  Gauss-Hermite  Quadrature  is  given 
in  Reference  [18]  as: 


E  = 


7 rm: 

*">771 


fKtuU't  2)  •  •  •  5  61 


(5.14) 


To  obtain  an  estimate  of  this  error,  the  divided  difference  (/  [G,  1 1,  ^2,  •  •  • ,  tm^mi  <fi] ) 

of  the  string  probability  function  is  obtained  by  evaluating  this  function  at  2m  +  1 
points  instead  of  m  points  as  required  by  the  Gauss-Hermite  Quadrature  rule.  This 
procedure  more  than  doubles  the  computation  time  required  to  calculate  the  proba¬ 
bility  of  a  given  string,  however,  this  is  necessary  to  control  the  error  and  arrive  at  a 
valid  decision.  Once  this  error  is  obtained,  it  is  compared  with  the  designated  error 
of  the  string  probability  function  and  the  total  error  is  taken  to  be  the  maximum 
of  the  two. 


String  Probability  Error  vs  Symbol  Probability  Error 

Both  algorithms  build  up  total  string  probabilities  from  individual  symbol  prob¬ 
abilities.  Since  the  noise  corruption  is  assumed  independent  between  measurements, 
any  given  string  probability  is  the  product  of  the  constituent  symbol  probabilities. 
Thus,  in  the  absence  of  estimation  errors  the  true  total  string  probability,  P(S),  is 
related  to  the  true  symbol  probabilities,  P(at)  by: 

;V 

P(S)  =  YlP(ai)  (5.15) 

;=i 

However,  exact  values  of  the  symbol  probabilities  are  not  available.  Thus  if 
we  define  the  total  string  absolute  error  to  be  E  and  the  ith  symbol  probability 
absolute  error  to  be  e;  then  equation  (5.15)  becomes: 

At 

P( S)  -f  E  —  JJ(P(a,)  +  e,)  (5.16) 

!=1 

As  stated  above,  total  string  probabilities  are  calculated  to  within  a  given  tolerance. 
In  order  to  guarantee  that  relative  error  on  the  total  string  probability  is  below  some 
maximum  (-Ej-clmax)’  a  maximum  relative  error  for  symbol  probability  estimation 
is  supplied  to  symbol  probability  estimation  routines.  The  defining  property  of  the 
relative  maximum  symbol  probability  error  is: 
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P{*i) 


—  erelmax’  ^  z  ^ 


(5.17) 


This  maximum  relative  error  for  the  symbol  probabilities  is  calculated  from  the 
designated  string  probability  relative  error  by: 


erelmax  —  V^relmax  +  * 


(5.18) 


This  guarantees  that  the  actual  total  string  relative  error  will  be  less  than  the 
designated  maximum  since: 
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(5.25) 


5.2.2  Conditional  Symbol  Probability  Calculations 

Estimation  of  the  conditional  symbol  probabilities  makes  use  of  standard  rou¬ 
tines  for  calculation  of  the  error  function  or  the  Marcum  Q  function.  Paris  method  [19] 
of  recursively  computing  the  Marcum  Q  function  was  used  for  single  level  crossing 
strings.  The  region  of  integration  for  the  octant  crossing  symbol  assignments  is  not 
the  standard  rectangular  region  usually  associated  with  a  bivariate  normal  random 
variable.  As  such,  approximations  to  the  region  are  made  by  sums  of  probabilities 
with  rectangular  regions  of  integration.  This  extension  of  the  bivariate  normal  cu¬ 
mulative  distribution  function  routine  is  considered  the  only  source  of  error  for  the 
octant  crossing  symbol  probability  routine. 


Octant  Crossing 


Octant  crossing  total  string  Cv..  -itional  probability  estimates  are  made  from  cor¬ 
responding  symbol  probability  estimates.  By  assumption  of  independent  Gaussian 
noise,  an  individual  octant  crossing  symbol  probability  given  a  particular  threshold, 
catalog  element  and  location  within  the  string,  may  be  estimated  by  the  formula: 


T  =  t,  Cj,  i) 


2na2 


'  )  dudv 


(5.2G) 


In  equation  (5.26),  12ai  is  the  region  of  plane  corresponding  to  symbol  a,.  The 
symbols  x,  and  y,  are  the  real  and  imaginary  component  of  the  ?'th  measurement. 

The  integral  in  equation  (5.26)  is  calculated  by  transforming  the  coordinates  so 
that  the  region  fiai  corresponds  to  the  7th  or  3rd  octant  (see  Figure  5.2).  Since 
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8 

270-360 

above 

Figure  5.2:  Octant  crossing  primitive  assignments. 

calculating  the  probability  of  a  rectangular  region,  of  a  Gaussian  bivariate  random 
variable  is  relatively  simple,  the  probability  of  the  7th  octant  is  approximated  with 
the  sum  of  probabilities  of  such  regions.  This  is  accomplished  by  first  determining 
a  lower  bound  for  the  integral. 

This  lower  bound  is  taken  to  be  the  integral  over  a  staircase  shaped  region  which 
approximately  covers  the  7th  octant  region  and  is  inscribed  by  the  7th  octant  region 
(see  Figure  5.3).  The  staircase  region  is  determined  such  that  the  points  of  the  stairs 
are  equally  spaced  in  angle.  The  probability  of  the  stancase  region  is  determined 
from  a  sum  of  the  probabilities  of  rectangular  regions.  An  upper  bound  for  the 
integral  is  then  calculated  by  summing  the  lower  bound  with  a  sum  of  probabilities 
of  rectangular  regions  between  the  inscribed  staircase  region  and  a  circumscribing 
staircase  region.  Thus,  the  resulting  region  integration  of  for  the  upper  bound 
probability  value  completely  covers  the  7th  octant  region  (see  Figure  5.4). 

If  the  difference  between  the  upper  and  lower  bounds  is  sufficiently  small  then  the 
total  7th  octant  probability  is  taken  to  be  the  average  of  the  two  values.  Otherwise 
the  angular  spacing  is  halved  and  the  binding  probability  measures  are  re-calculated. 
The  probability  of  the  3rd  octant  is  calculated  in  a  similar  manner. 
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Figure  5.3:  Inscribed  staircase  region,  5  points 


Figure  5.4:  Circumscribed  staircase  region,  5  points 
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Single  Level  Crossing  Conditional  Symbol  Probability  Calculation 


For  a  fixed  threshold,  the  definition  of  the  single  level  crossing  scheme  indicates 
that  the  probability  of  an  entire  single  level  crossing  string  may  also  be  calcu¬ 
lated  as  the  product  of  probabilities  of  a  series  of  Rician  random  variables  being 
above  or  below  the  given  threshold.  For  simplification  of  the  program,  single  level 
crossing  strings  are  transformed  into  a  strings  of  binary  symbols  in  which  each  sym¬ 
bol  indicates  “above”  (1)  or  “below”  (-1).  These  symbols  are  thus  considered  the 
constituent  symbols  of  the  strings  and  total  string  probability  (given  the  catalog 
element  and  the  threshold)  in  equation  (5.2),  etc.  is  estimated  by  the  product  of  the 
probabilities  of  a  series  of  Rician  random  'ariables,  all  exceeding  or  falling  below 
the  given  threshold.  Thus,  the  probability  of  a  transformed  single  level  crossing 
symbol  is  “above”  the  threshold,  given  a  particular  threshold,  catalog  element  and 
location  within  the  string  may  be  estimated  by  the  formulas: 


P(a,  |  T  =  t,Cj,i)  = 


(5.27) 

(5.28) 


In  equation  (5.28)  and  <72  are  as  defined  above  and  Qj  is  the  Marcum  Q  function. 

The  “below”  case  is  calculated  using  different  integration  limits  or  the  comple¬ 
mentary  Marcum  Q  function.  Total  string  probability  is  calculated  from  products 
of  these  symbol  probabilities. 


5.2.3  Conclusions  and  Results 

Comparison  of  likelihood  functions  generated  with  the  analytical  methods  and 
those  generated  by  the  histogram  method  are  similar  for  the  noise  levels  which 
have  been  run.  Misclassification  results  are  given  in  Figures  5.5  and  5.6.  These 
figures  represent  results  obtained  under  the  assumption  of  completely  known  noise 
variance.  It  is  expected  that  the  same  matched  noise  phenomenon  exhibited  by  the 
single  level  crossing  scheme  under  the  histogram  method  will  appear  if  mismatched 
noise  levels  are  used. 
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Figure  5.5:  Octant  crossing  vs.  coherent  nearest  neighbor 

For  comparison  purposes,  results  from  competing  algorithms  which  employ  the 
entire  measurement  vector.  While  the  performance  of  the  nearest-neighbor  classi¬ 
fiers  is  usually  better  than  that  of  the  level-crossing  classifiers,  the  level-crossing 
based  pattern  representations  can  be  specified  with  only  ^  to  the  number  of 
bits  required  by  the  nearest -neighbor  specification.  Thus,  while  the  level-crossing 
observations  cannot  be  considered  “sufficient  statistics’’  they  are  highly  information 
preserving. 

The  consequences  of  this  new  algorithm  are  far  reaching.  Inference  of  more 
complex  generating  grammars  requires  the  use  of  an  informant  [20].  The  informant 
in  these  algorithms  is  expected  to  answer  queries  about  the  applicability  of  any  given 
string.  With  the  histogram  method  of  likelihood  function  estimation,  knowledge  of 
the  category  of  a  given  string  is  dependent  upon  the  conditional  probability  of 
occurrence  of  that  string.  If,  for  all  catalog  elements,  a  string  occurs  with  very 
small  probability  then  it  may  take  a  very  long  time  for  knowledge  of  the  category 
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Figure  5.6:  Single  level  crossing  vs.  nearest  neighbor 

to  become  available.  The  fact  that  any  arbitrary  string  may  be  classified,  no  matter 
how  small  the  conditional  probability  of  occurrence,  means  that  automatic  inference 
of  more  complex  generating  grammars  is  possible  with  this  technique. 

The  ability  to  classify  a  given  symbolic  representation  quickly  and  easily  also 
provides  more  flexibility  in  terms  of  classification  algorithm  design.  The  major  hin¬ 
drance  in  the  determination  of  the  symbolic  representations  conditional  probability 
is  the  dependence  on  the  adaptive  magnitude  threshold.  It  is  proposed  that  the  in¬ 
vestigation  of  non-parametric  syntactic  pattern  recognition  techniques  include  the 
investigation  of  symbolic  representations  which  do  not  include  adaptive  thresholds 
which  are  calculated  over  the  entire  range  of  measurements.  Such  symbolic  assign¬ 
ment  algorithms  could  represent  quantities  such  as  differential  phase  or  magnitude. 
Furthermore,  inference  of  the  generating  grammar  for  such  symbolic  assignments 
could  be  made  directly  from  noiseless  measurements. 

It  is  conceded  here  that  there  exists  no  reason  to  believe  that  the  generation  of 


91 


such  ad-hoc  symbolic  representations  can  be  accurately  modeled,  in  general,  with  a 
given  generating  grammar.  However,  by  restricting  the  range  of  measurements  and 
using  automatic  methods  of  inference,  a  simple  classifier  which  uses  an  extremely 
distilled  statistic  of  the  full  observation  has  been  demonstrated  to  give  good  per¬ 
formance,  even  when  compared  to  a  classifier  using  the  original  vector  observation. 
Continued  exploration  in  this  area  may  provide  insight  into  the  “structure”  of  such 
measurements. 

As  a  further  note,  computation  time  for  these  strings  can  be  exceedingly  long. 
This  becomes  true  when  string  probabilities  are  close  for  two  or  more  catalog  ele¬ 
ments.  Much  effort  has  been  expended  on  increasing  computational  efficiency.  Some 
improvements  have  come  about  from  novel  storage  and  retrieval  schemes  of  string 
and  individual  symbol  probabilities. 

5.3  Conversion  of  Syntactic  Classifiers  to  Deterministic  Forms 

The  syntactic  classifiers  derived  in  [17]  are  non-deterministic  in  nature.  This 
means  that  there  may  be  more  that  one  production  rule  which  is  applicable  to  a 
particular  nonterminal  in  a  sentential  form.  In  terms  of  automata,  the  transition 
function  may  have  more  than  one  value  for  a  given  state  and  symbol  (the  automata 
may  have  more  than  one  place  to  “go”  on  a  given  symbol).  Deterministic  forms  of 
automata  or  grammars  do  not. 

Derivation  of  a  deterministic  form  for  the  syntactic  classifier  from  the  non- 
deterministic  form  can  be  summarized  as: 

1.  Label  the  stop  state  each  catalog  elements  automata  with  the  corresponding 
catalog  element  name. 

2.  Merge  the  non-deterministic  automata  with  labeled  stop  states  together  by 
adding  a  new  start  state  with  e  (null  string)  transitions  to  the  start  states  of 
each  of  the  automata. 

3.  Convert  the  resulting  non-deterministic,  merged  automata  to  a  deterministic 
form. 
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The  conversion  to  a  deterministic  form  is  a  well  known  process.  It  is  fully 
describe  in  [21].  One  of  the  more  noteworthy  aspects  of  the  process  is  that  each 
state  of  the  deterministic  form  of  the  automaton  is  equivalent  to  a  set  of  states  of 
the  non-deterministic  form  of  the  automaton.  This  leads  to  the  fact  that  there  are 
potentially  2n  states  in  the  deterministic  form  of  the  automaton,  where  n  is  the 
number  of  states  in  the  equivalent  non-deterministic  automaton. 

Furthermore,  the  stop  states  in  the  deterministic  form  are  determined  from  the 
stop  states  of  the  non-deterministic  form.  Specifically  a  state  in  the  deterministic 
form  or  the  automaton  is  a  stop  state  for  catalog  element  i  if  th°  stop  state  for 
catalog  element  i  is  contained  in  this  state.  Since  states  of  the  deterministic  form 
of  the  automaton  correspond  to  sets  of  states  from  the  non-deterministic  form  it  is 
easy  to  see  that  stop  states  can  potentially  correspond  to  acceptance  of  a  string  by 
more  than  one  catalog  element.  This  is  the  result  which  allows  for  analysis  of  the 
grammatical  inference  procedure. 

Even  though  application  of  the  algorithm  for  derivation  of  the  deterministic 
automata  has  potential  to  increase  the  number  of  states  exponentially,  it  has,  as 
a  rule,  decreased  the  number  of  states.  Application  of  the  algorithm  to  derived 
automata  has  revealed  a  number  of  multiply  accepting  states  are  present  in  the 
classifier.  These  states  may  be  “traced  back”  to  determine  exactly  which  strings 
are  being  accepted  in  a  multiple  sense. 
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(a) 


AJ.C  C 


(b) 


Figure  5.7:  Conversion  to  deterministic  form,  (a)  non-deterministic  automata  set 
with  labeled  final  states,  (b)  merged  non-deterministic  automata,  (c)  resulting 
deterministic  automata 
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